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ABSTRACT 

We investigate the role and behaviour of dust grains in C-type MHD shock waves 
in weakly ionized, dense molecular clouds. The structure of C-type shocks is largely 
determined by the coupling of the charged species and the magnetic field. In weakly 
ionized clouds, charged dust grains enhance the energy and momentum transfer be- 
tween the magnetic field and the neutral fluid, and dominate the neutral collisional 
heating rate. New shock models are developed here for steady oblique C-type shock 
structures with shock speed v s — 18km/s, pre-shock number density n H = 10 5 cm~ 3 , 
and a grain population represented by either a single grain species or a MRN grain size 
distribution. The grain size distribution is calculated using Gauss-Legendre weights 
and the integrals over the continuous distribution of grain sizes are converted to a 
series of separate grain bins or 'size classes'. The dynamics of each grain size class 
various through the shock front; smaller grains remain coupled to the magnetic field 
and larger grains are partially decoupled from the magnetic field due to collisions with 
the neutrals. The charges on the grains are allowed to vary, via the sticking and re- 
releasing of electrons, increasing with increasing electron temperature. The increase in 
grain charge increases the coupling of the grains to the magnetic field, and magnetic 
field rotation out of the shock plane is suppressed. MRN(mantles) and MRN(PAHs) 
distributions are also compared with the standard MRN model. Increasing the grain 
sizes in the MRN(mantles) model leads to an increase in the collisional heating of the 
neutrals leading to hotter, thinner shock structures than those using a standard MRN 
distribution. With the addition of PAHs, the electron abundance is reduced and the 
grain charge is held constant, resulting in less grain coupling to the magnetic field, and 
substantial rotation of the magnetic field out of the shock plane. The effects of varying 
the orientation of the pre-shock magnetic field Bo with the shock normal, specified by 
the angle 9, are also considered. It is found that there are critical values of 8 below 
which the shock is no longer C-type and the transition becomes C* or J-type. The 
degree of non-coplanarity of the shock solution depends upon the grain model chosen, 
as well as the angle 6. 
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1 INTRODUCTION 

Interstellar shock waves play an integral part in the chemical 
and physical evolution of the interstellar medium. They are 
driven by cloud -cloud collisions, stellar winds and supernova 
explosions fe.g. lMcCrav fc Snowll979T) . They cause large in- 
creases in the local pressure leading to both compression and 
heating of the gas. As a result, the local increase in density 



* E-mail: jchapman@physics.mq.edu.au 
f E-mail: wardle@physics.mq.cdu.au 



leads to an increase in the rates at which chemical reactions 
can occur. C-type sh ocks may be associa ted with strong in- 
frared emission lines iDraine et alJll983ft . and shock models 
have been used to explain li nes in the K-L region of Orion 
JPraine fc Robergd Il982t IChernoff. Mckee fc Hollenbachl 
Il982fl . C-tvpe shocks may also be responsible for the ab- 
sorbtion lines along lines of sight, in particular of CH+, in 
diffuse clouds t Flowcr.Pinemi^e^F^oret^^Hartauist Jll985l : 
iDraine fc Kat dll986al lb UPineau des Forets et alJll986lh 

The nature and structure of shock waves travelling 
through molecular clouds are strongly dependent upon the 
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strength of the ambient magnetic field Bo as well as the local 
ionization fraction. When the magnetic field strength is large 
and fractional ionization low, a multifluid magnetohydrody- 
namic (MHD) description of the fluid is usually required 
<MuUarJll97l . In the multifluid treatment, there are two 
speeds at which compressive disturbances may propagate. 
In the charged component of the medium, disturbances can 
propagate at the fast magnetosonic speed which is depen- 
dent upon both the magnitude |_Bo| and the orientation of 
Bo with the local shock normal. In the neutral fluid, com- 
pressive disturbances propagate at the neutral sound speed. 

The charged fluids gain energy and momentum from 
the magnetic field B, and are coupled to the neutral fluid 
via collisions, transferring energy and momentum, acting to 
weakly couple the neutrals to B. If the fractional ionization 
rate is small enough, leading to a low neutral-ion collision 
rate, the two fluids c an be thought of as separate an d in- 
terpenetrating fluids llDraindll98(il : iDraine et alJll983ft . The 
ions and electrons contribute only small amounts of mass 
and pressure, however the magnetic fields are strong enough 
to influence the shock dynamics. If the shock speed is less 
then the (fast) magnetosonic speeds of the charged fluids, 
the charged species are accelerated ahead of t he neutrals, vi a 
ambipolar diffusion in a magnetic precursor jDraindll98Cl) . 
If the neutral temperature remains low enough inside the 
shock and the neutral flow is everywhere supersonic, the 
shock is C-type and all the hydrodynamic variables remain 
continuous. 

At low enough densities, a three fluid system is ad- 
equate with the inclusion of neutrals, ions and electrons. 
The grains may be ignored since their contribution to the 
energy and momentum of the neutrals is negligible. For 
densities iih ~ 10 2 cm -3 the energy and momentum trans- 
ferred from the charged d ust grains to the ne utrals can no 
longer be neglected (e.g., IDraine et al.lll983l) . and in re- 
gions of high density, riu ~ 10 5 cm -3 , and the grains domi- 
nate the collisional heating of the neutrals. The grain abun- 
danc es in the interstellar medium obey a size distribution 
fe.g..lMathis. Rumpl fc Nordsiecrfl977tlDraine fc Leell984t 
IWeingartner fc Draind 1200 ill , and for large grain sizes the 
friction between the grains and neutrals can become great 
enough that the grains are decoupled from B. In the oppo- 
site limit, smaller grains may remain tied to B. Dust grains 
also have the potential to drift with respect to the neutral 
fluid and move out of the shock plane, since grain-neutral 
collisions can cause their partial decoupling from B leading 
to non-coplanar shock structures. 

For perpendicular shocks with Bo oriented perpendic- 
ular to the local shock normal n, changes in B and the 
dynamic variables inside the shock ar e restricted to the 
shock plane containing Bp and n (e.g., IDraine et al.lll983 
iFlower fc Pineau des Fore"tsl [2003) . When Bo is oriented 
obliquely with the shock normal, there is a current paral- 
lel to the B component transverse to the shock propagation 
direction. Bi, the B component perpendicular to n, can 
rotate out of the shock plane inside the shock before return - 
ing to the shock plane downstream. IWardle fc Draind (Il987t) 
demonstrated that by changing only the transverse compo- 
nent of B, very different hydrodynamic shock structures are 
found. 

iPilipp fc Harauistl l|l994fl investigated the rotation of B 
in fast steady shock structures with the inclusion of one grain 



species. This involved integration of the shock fluid conserva- 
tion equations from upstream to downst ream which proved 
difficu lt for shock speeds ~ Skms" 1 . IPilipp fc Harauistl 
( 1994) concluded that st eady fast shoc k solutions do not 
exist for these conditions. IWardld il998fl showed that these 
solutions were in fact intermediate, not fast, solutions which 
do not exist for shock speeds above a critical value. IWardld 
(1998) presented magnetic field (B x versus B y ) phase space 
topologies, demonstrating that the shock solution belongs to 
a one-parameter family. Integration from upstream becomes 
complicated since the phase space trajectories can diverge 
due to finite numerical precision, and end up on neighbour- 
ing intermediate solutions which can become unphysical. In- 
tegrating from the downstream state ensures that the correct 
fast solution trajectory is followed. 

In this paper we investigate the non-coplanar nature 
of C-type shocks with the inclusion of a grain size distri- 
bution. The shock is assumed steady state since the time 
scales over which the fluid variables vary are longer then the 
time required to transverse the shock. The methodology of 
Wardle (1998) is extended to include a grain size distribu- 
tion and the integration of the energy equation allows for 
radiative cooling. The shock problem comprises of three or- 
dinary differential equations (ODEs) in B x , B y , and pressure 
P, along with algebraic expressions for the neutral velocity 
components, electric field components, and charged species 
properties. It is a two point boundary problem connecting 
the upstream and downstream states, and integration pro- 
ceeds from downstream to upstream to ensure the fast shock 
solution. With radiative cooling the fluids do not cool to 
their final temperature until far downstream of the shock 
front , and in tegration over such long time scales is infeasible 
fe.g.. lDrainelll98Cl) . A shooting integration method is imple- 
mented which searches for an initial state inside the cooling 
zone which, after integration towards upstream, yields the 
pre-shock state. Conditions inside molecular clouds are as- 
sumed, uh = 10 5 cm -3 , with shock speed v s ~ 18 km s _1 . 

The derivation of the governing shock equations are 
given in Section |21 the shock jump conditions are derived 
in Section |3J and the heating and cooling processes are dis- 
cussed in Section^] The treatment of the charged species are 
addressed in Section |SJ along with the derivation for both 
the electron and grain Hall parameters and the electric field 
component E' z , in Sections |S| and |7| respectively. The shock 
calculation and integration methods are explained in Sec- 
tion |S] followed by the results in Section |U] A discussion and 
summary are presented in Sections 1 1 01 and [Til 



2 FLUID EQUATIONS 

The fluid is weakly ionized, and the inertia and thermal pres- 
sures of the charged species are neglected, as well as most 
processes that create and destroy species. There is only the 
one exception, when grain charging is used (equation 15611 ) 
the electrons are allowed to stick to the grains, and so elec- 
trons are removed from, or injected into, the electron fluid 
to conserve charge density (equation 12411 ). Each species is 
characterized by a mass in, mass and number densities p and 
n, respectively, charge Ze and velocity v. Charged species 
are given a subscript of either e, i or g, denoting electrons, 
ions and grains, respectively. A subscript of j represents 
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any charged species. The pre-shock medium has mass and 
number densities po and no, respectively, magnetic field Bo, 
pressure Po, temperature To and the incoming fluid has a 
speed v s . The shock is assumed steady and plane-parallel 
d/dt = d/dx — d/dy — and the coordinate system is de- 
fined in the shock frame with Bo lying in the x — z plane 
making an angle 9 with the z axis. The incoming flow is 
parallel to the z axis. The equations for a steady multifluid 
flow are then given by the conservation of mass 



The current density is given by 



dpv z 



= 0, 



J x B 



dz 

momentum 

dv dP 
PVz-r- + -7-z = 

dz dz 

and energy 

d(Uv z ) pdv z 
dz dz 

along with the following from Maxwell's equations: 
dEx . v q 



G- A, 



dz 
and 
dB z 
dz 

Lastly, from Ampere's law: 



0. 



dB. 

dz 
dB y 

dz 
and 

J z = 0. 



47T 

Jy-, 

dz c 

4n 

O X 

C 



(i) 

(2) 

(3) 

(4) 

(5) 

(6) 
(7) 

(8) 



P and U denote the neutral thermal pressure and internal 
energy, G and A, the heating and cooling rates per unit 
volume respectively, and E, B and J are the electric field, 
magnetic field, and current density, respectively, with 

pkT 



P = 



(9) 



The internal energy U is the sum of the thermal kinetic 
energy and excitation energy of the excited molecules 



(10) 



where the summation is over the electronic, rotational and 
vibrational levels I, with corresponding probabilities pi and 
energy Ei. In a H2 gas the vibrationally excited levels have 
negligible populations, and the rotational levels have a ther- 
malized distribution, and the sum can be taken as unity 
JWardld Il99ll: lHartauist fc Casellil Il998ft which is a good 
approximation in the hottest part of the shock where the 
pressure becomes important. Equation 1101 becomes 



U = 



7-1 

where 7 is the adiabatic index with 7 = 7/5. 

From mass conservation of the charged species 



d 

Tz^ 1 



1) = 0. 



(ii) 



(12) 



1 ^3 ( v j 



assuming charge neutrality 



(13) 



(14) 



(vj — v) is the drift velocity of a species j with respect to 
the neutrals. The grains obey a continuous size distribu- 
tion implying equation 1141 contains an integral over grain 
size. However, the MRN size distribution is calculated using 
Gauss-Legendre weights (Section 15. 2^ . and the distribution 
is represented by a series of discrete grain size bins each 
with a specified number density and charge. Each grain bin 
is then considered a separate species j in equation I)1M|I . 

The rotation of B±(|B±| = y/B% + Bf) inside the 
shock front is driven by the Hall current along the E'xB 
direction, and the Pedersen current along the E' direction. 
The difference in drift speeds between charged species drives 
the Hall current. In the ambipolar diffusion limit all charged 
species are highly coupled to B and the Hall current van- 
ishes. In this limit the magnetic forces on the charged species 
dominates the neutral collisional drag which is the case for 
the electrons and ions in molecular clouds. When charged 
dust grains are present the neutral-grain drag force is non- 
negligible resul ting in a non-zero H all current and there is 
rotation of B 1 . IWardle fc Nel l)l999f) demonstrated that for 
dense molecular gas, the contribution of dust grains to the 
Hall conductivity leads to substantial changes in the dynam- 
ics of the gas as well as the evolution of B . 

Treating the charged particles as test particles, and ne- 
glecting their inertia, their drift through the neutrals is 



E+^xB 



) +7jPjP(v-Vj) = 0, 



(15) 



which represents the balance between the electroma gnetic 
forces and the collisional drag ijDraine! ll98rtlShJl9"83T) . The 

rate coefficient for elastic scattering between particles of a 
species j and the neutrals, < av >_, , is related to jj via 



li = 



rrij + m 



(16) 



If a is not a function of v, < av >j may be written as (T jUj, 
where Uj is the effective velocity given by llDrainelll98ri 



= [Vi + 



where 



<Pi = 



128 ( kT kT. 



97T 



+ 

mjv m 



(17) 



(18) 



The treatment of < av >j for each charges species is dis- 
cussed further below in Section 

Equation 11511 can be recast in terms of the electric field 
E' in the frame co-moving with the neutrals, 



where 



(B ■ E')B 0j E'xB 



B'< 



+ 



B 2 

B x (E' x B) 



S3 



E' = E + v x B/c 



(19) 



(20) 
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and f3j, the Hall parameter for species j, 

a _ ZjeB 1 
Pi 



rrijC jjp 



(21) 



is the product of the gyrofrequency and the time-scale for 
momentum exchange with the neutral fluid. 

Integrating the mass conservation equations Q and 
11211 along with charge neutrality 1141 1 gives (the subscript k 
runs over the ions and grains): 



pv z = pov s 

pkV z k = POkVs 

n e = m + S~] Zg 



(22) 
(23) 
(24) 



where a local equilibrium approximation is used for Z g 
(equation 15611 . Combining equations with equations 
JSJ, and integrating, gives the neutral velocity components: 



and 



1 

Pov s 
1 

Pov a 
1 

pOVs 



4ti 



(B x — Box) 



4ti 



P -P + p v 2 + —(B% x 

07T 



B x — B, 



(25) 
(26) 

(27) 



Equations @ and © are combined, along with <f 1 1 p . to 
construct an ODE in P: 



dP 

dz 



1-% 



( 7 -l)(G-A) 



+ 



7 P 



a x — h ±f „ — — 

da az 



where c s is the neutral sound speed 

Equation 1281 breaks down at the sonic point v z 
Finally, from equations Q and (JSJ, 



-Boa: 



E x =0, 

c 

and 

S z = Bo z — cos 

Integrating equation 1201 . and using 13011 : 



Si = (uj/Bz - 
and 

Ey = (— VsBqz 



V Z By 



v,B x - v x B, 



1 



(28) 
(29) 

(30) 
(31) 
(32) 
(33) 



Charge neutrality (equation 1141 is used to find E' z (see Sec- 
tion 0. 

The shock problem then comprises of three ODE's in 
B x , By, and P equations ©, 0, and 1281 1 . respectively, 
with the algebraic relations 1221 . I251 - I27|l for the neutral 
density and neutral velocity components, respectively, along 
with equations l31H - 133t for B z , E' x and E' y . To calculate the 



right hand side of the ODE's © and (J7J, the current den- 
sity components J x and J y need to be known, which are 
dependent upon the drift velocities, Vj — v (equation 1131 1. 
Calculating Vj — v from equation 1191 requires knowledge of 
the Hall parameters j3j, as well as the electric field compo- 
nent E' z . The treatment of f3j and E' z are described below in 
Sections [S] and [7] The heating and cooling rates needed 
for the energy equation 1281 are discussed in Section [1] 



3 JUMP CONDITIONS 

The boundary conditions for equations @ , Q , and 1281 are 
that the derivatives vanish far upstream and downstream of 
the shock. Shock solutions begin and end at points where 

dB x _ dB y _ dP _ Q 

dz dz dz 



(34) 



From equations JSJ-©) J = 0, and equations 1 1 ■'->> and l|191 

give E' = and Vj — v = 0. All species are co-moving 
upstream and downstream of the shock and obey the same 
overall jump conditions. Also Equation 1281 gives G = A. 
E' = 0, so E = —v x B/c and by equations i3"2l and ll3^l : 



V Z By = 



(35) 



and 



v z B x - v x B z = v s B 0x . (36) 
Eliminating v y from equations 1261 and 1361 gives either 



By 



or 

Bl 



(37) 
(38) 

47rpolIs 

The second solution implies v z is equal to the intermediate 
speed, resulting in a rota tional disconti nuity and no com- 
pression across the shock jCowlindll97rJ) . So taking B y = 
gives the condition v y = ahead of and behind the shock. 

In a radiative shock, the downstream fluid radiates away 
the energy transferred to it by the processes inside the shock. 
Cooling proceeds at a rate A(p, T) (Section ^} over a post 
shock cooling layer. Drift of the charged particles inside the 
shock lead to increases in thermal energy due to elastic scat- 
tering between the charged and neutral components. The 
immediate post-shock temperature is therefore much higher 
than the final post-shock temperature after the fluid has ra- 
diated away the energy inside the post-shock cooling layer. 

The fluid cools inside the cooling layer at an almost 
constant pressure as the density increases and fluid speed 
decreases. The post-shock fluid cools until it returns to the 
pre-shock temperature To. Applying the isothermal jump 
conditions T = To along with equations JUJ and 1221 gives 

£ = - = ^. (39) 
Pa n v z 

Equations 1251 , 1271 , 1361 , and 1391 give a cubic in v z / v a 



+ ao-| + 01 h a 2 = 0, 



where 



ai = 



M 2 A0 



(3p 2 
2u 2 



(1 



JM 2 



M 2 A0 



(40) 

(41) 
(42) 
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and 



«2 



(43) 



Only roots with < v z /v s < 1 are suitable for there to be 
compression across the shock. Once the downstream v zc i is 
calculated, the corresponding P, B x , and v x are easily found. 
One could alternatively obtain a cubic in B x and firstly solve 
for the downstream component B xd instead. 



4 HEATING AND COOLING PROCESSES 

The main heating processes inside dense molecular clouds 
are the heating by cosmic ray s and heating from shocks . 
The cosmic ray heating rate is iGoldsmith fc Langej|l97^) 



G cr = 6.4 x 10 n(H 2 ) er S s cm s " 



(44) 



Elastic scattering between two species at different tempera- 
tures will lead to heat exchange, as well as heat generation 
if the species are streaming relative to each other. The heat- 
ing rate per unit volume prod uced by elast ic scattering of a 
species a with a species ft, is foraindll986h 

G ali = , Papp , 7 a g[3k(7> - T a ) + m \v a - v„| a ], (45) 



(m a + mg) 



where 



7«/3 



< OV >a/3 



(46) 



If a species j is unable to cool efficiently and has a small heat 
capacity per unit volume then G jn = (|Chernofl1ll98if l and 



(47) 



T i =T+—m n \v j 
The heating rate of the neutrals due to species j is then 
G" i =Pft7j|v-v J | a . (48) 
Thus, by equation 1151 . the total heating rate is 

g " = pE^! v - v ;I 2=j ' e '' ( 49 ) 

3 

The low mass of the electrons ensures that their collisional 
heating is negligible. With low fracti onal ionisation th e 
grains dominate the collisional heating jDraine et aljfl98^l . 
The total neutral heating rate is then G = G™ + G cr - The 
dominant cooling mechanisms inside shock-heated regions 
is collisional excitation of molecular H2. Neglecting the vi- 
brational cooling, which is negligible for ri(Ha) ~ 10 7 cm -3 
T ~ 3000K, we adopt the cooling rate 

L r HLrL 



A = n(H 2 



+ L r 



ergs cm 3 s , 



(50) 



with the rotational cooling rate coefficients for high and 
lo w density, L r H and L r L, from equations (10) and (11) 
of lLepp fc Shu3<1983ft . 



has a 1/v dependence for drift speeds below 20km/s, and 
the rate coefficient for ion-neutral scattering is < at) >;~ 
1.6 x 10 _9 cm 3 s _1 . The ion Hall parameter (equation 1)2111 ) 
can thus be written in terms of the pre-shock Hall parameter 



R - R B Vz 
Pi — PiOS • 

Inside a molecular cloud (lElmegreenlll979l) 

Tin V 10"cm" 



(51) 



(52) 



The ions are unable to cool efficiently, as they have a small 
heat capacity per unit volume, due to the low fractional 
ionisation, thus G nl w and equation (1471 applies to T;. 
The electron scattering cross-section is a e ~ 1 x 10 _15 cm 2 , 
and as m e << m n , then by equation 1171 : 



U e = Up e + |v e — V 



21 2 



where 



fe = 



128 kT e 
9ir m e 



(53) 



(54) 



The electron temperature T e is needed to calculate tt e 
and Z g . Electrons and ions are highly coupled to B with 
\fti, e \ » 1 (equation 1211 1. and thus (vj — v) w (v e — v) w 
E' xB/B 2 (equation ([T^ll. Equation lETTl. therefore, implies 
that T e ~ Ti, however IPraine et al.l dl^831 showed that in- 
side shocks in dense molecular gas the electrons lose energy 
due to impact excitation of H2 and T* e <Tj. T e is t herefore 
approximated via (see Fig. 1 of lDraine et alJ il983f) 1 



T e «T + 0.2(Ti -T). 



(55) 



Once T e , and subsequently u e , are known, then ft e can be 
determined. The calculation of ft e is complicated by the de- 
pendence of u e on |v e — v|. To determine ft e , equations 1171 . 
1191 . and 1211 must be solved, this is discussed in Section|S| 



5.2 Grains 

The charge on a dust grain Z g results from collisions with 
other charged particles as well as the ejection of photoelec- 
trons. The charge state of a grain contains stochastic fluc- 
tuations about some mean JSpitzerlll97l iGail fc Sedlmavrl 
Il975l) . If the fluctuations are neglected, the instantaneo us 
mean charge on a dust grain of radius a g is jDrainelll980T) 



(56) 



Equation l|56|l is only valid for \Z g \ S> 1 which is not al- 
ways satisfied for the pre-shock conditions of interest; in 
cold molecular clouds T e o = 10K and for a MRN distribu- 
tion a g = [50A, 2500A], \Z g0 \ < 1. We force Z g = -I when 
\Z g \ < 1 by equation 1561 . The scattering cross-section for a 



spherical grain is a g 



and by equation 1171 , 



5 CHARGED SPECIES 
5.1 Ions and Electrons 



The ions are sing l y cha rged with mi = 30mH- As dis- 
cussed in IWardld (1998), the scattering cross-section <t; 



128 
~97r~ 



kT 
rriN 



(57) 



as the grain thermal velocity is negligible compared with 
that of the neutrals. The internal density is 2.5g/cm 3 . By 
equation 1561 the relative gyrofrequency (Z g eB /m g c) for a 
smaller grain a g3 is greater than that of a larger grain a g i 
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by (a g i/a gs ) 2 , and so the smaller grains remain better cou- 
pled to the field. Consequently, the smallest grains have the 
largest f} 9 . The calculation of g is addressed in Section 

The abundance of the pr e-shock grains is given by a 
MRN grain size distribution llMathis. Rumpl fc Nordsiecld 
Il977fl which assumes a mixture of silicate and graphite; 



dn g o , -3.5 
— — = An H a , 
da 



A = 1.5 x 10" 25 cm 2 5 . 



(58) 



The limits are [m, a z ] = [50A, 2500A] and defining y = log a; 

floga 2 



UgO 



Aa 3 ' 5 da 



Aa °dy. 



(59) 



log ai 



Inside dense molecular clouds, the grain size distribu- 
tion may differ cons iderably from the MRN model (e.g. 
iLi fc Greenber3l2003D . Observational measurements of the 
visible extinction of objects behind dense clouds su ggest 
that larger dust grains may be present (e.g. IWhittet et ail 
11983ft . Inside dense clouds, the growth of ice mantles and 
dust coagulation can remove the smaller grains (a g « 
IOOO A), and can increase the grain mass by a factor of up 
to 2 JChokshi et al]ll993l) . We also consider a MRN dis- 
tribution with increased limits [aim , 02ml = [90A , 4500 A] 
( Nish i~Nakano fc Umebavashill99 if) . The distribution is re- 
normalised with y4 raa „nesflm ' 5 = Aa~ 2 5 so A mantle3 = 
A(a/a m )- 2 - 5 =6 ,^'^- 25 "~ 2 - 5 
sity is 1.14 g/cm 3 

Using Gauss-Legendre weights (a highly accurate 
method of integrating smooth functions) the MRN grain 
size distribution may be represented by the summation 
over a number of discrete grain bins or grain size classes 
iPress et al.lll993) and equation I59H may be expressed as 

/■logai N N 

xao = / f{y)dy = 53 Wmfiym) = 53 x g0m , (60) 

" log a 2 



6.5 x 10 cm ' and the mean interior den- 

3 



m — 1 



m— 1 



-2.5 



where y m are abscissae with weights w m , f{y) = Aa~ 
and y — log a. The y m and w m are pre-determined from 
the Legendre polynomial of order N. Each of the N grain 
bins, have an associated size a, number density x g o m , and 
corresponding parameters, Z gm p gm , v gxm , v gym , and v gzm . 
Using equation 16011 . any integrals over the grain size distri- 
bution can be generalised such that for any function g(y); 



log ax 



then the drift velocity equation 1191 may be written as 
P 2 E' ± x B p. 



v j - v _ Rv > 
— ~ — - "i^W 



l + P 2 B 2 + 1 + P. 



rE 



(65) 



Taking the dot product of Vj — v with itself gives 

I |2 o2m'2 2 | Pj t-,'2 2 

|vj-v| =P j E ][ c + T -^-E x c. 

Defining otj as 

ZjeB m, + m 

&j = , 

rrijc (jjp 

then equation 12 II becomes 

= Z je B (mj+rn) = Oj_ 
J rrijC GjUjp a, ' 

Eliminating Uj from equation 168H using 1171 gives 



ft 



- Pi- 



rn 



(67) 



(68) 



(69) 



After eliminating |vj — v| from equations 1661 and 169H then 

^ E ;, 2 + ^(E; 2 + E; 2 + f) 

(70) 



+ P M.2 ,2 „2 U ' 



which has only one positive root for /3 2 . The root — \Pj\, is 
then taken for the negatively charged electrons and grains. 



7 CALCULATING E' z 

Lastly, E' z is needed, and using equa tions (|14| an d I23H . 
equation l|19fl is recast in terms of Vj z l|Wardlelll99Sft : 



v jz = n j0 Z j (p j E' z + qj), 
where 



log a 2 



f{y)g{y)dy= 53 w mfijjm)g{y m ) = 53 x amg{y m ). (61) 



m — 1 m — 1 

Equations 1491 and 1131 1 are calculated in this way, for ex- 
ample; 



Pi = 
and 

1i = 
+ 



c p 3 (B 2 p 2 +B 2 



n j0 Zj 1 + p 2 



rijoZj 
1 



B J - 



cp 2 (E' x B) 2 



(71) 



(72) 



1 + P 2 B 2 



n j0 Zjl + p 2 



(B X E' X + ByEy)-^. 



J = eriH 



s v e + 53 

m — 1 



(62) 



From equations JTIJ, and JTTt 

n e = ^Z k n fc — = ^-_ 



Qk 



(73) 



(74) 



6 HALL PARAMETERS 

To determine p e and P g , equations I17H . 119L and 12 II are 
solved as follows. Writing 



e; = 



and 



E', = 



B 3 



B x (E' x B) _ E' (B ■ E') 



B 3 



B 



B'< 



B, 



(63) 
(64) 



where k runs over the grain and ion species. Also by equa- 
tions © and (EU; 



= — n e v ez + 53 Z k n k 



v z k = —n, 



J2z k n k0 . (75) 



Using equation I71H for the electrons with equation 17511 . 
and using equation l|740 to eliminate n e gives 



_1 n k pZ k 1 



(p e E' z + q e ) A r d rieo 



p k E' z + q k 



= 0. 



(76) 
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tions for E' z . In lWardld(ll99gl) . N = 3 and equation C5) is an 
easily solvable quadratic in E' z . For arbitrary N the choice of 
the correct root for E' z is not so simple. Imposing the require- 
ment E' z = upstream and downstream enables the identi- 
fication of the bracketing pole species (—qjo/Pjo)- Since E z 
is continuous through the shock, the correct E' z root will 
remain between these bracketing species poles. Since 



PjO 



n j0 Z j0 l+Pj 



1 



So 



R2 
D 



rijoZjo 



(77) 



(78) 



then pjo < and the sign of qjo is determined by the sign of 
Zjo- The ion pole —qio/pio will always be negative, however 
poles of the negatively charged species are positive, thus 
—qio/pw will always be the minimum bracketing pole. /3 e o ~ 
10 3 and f3 g0 ~ 10 1 - 10 2 for the grain sizes considered, so 
\qeo/peo\ < \q g o/p g o\- Thus the solution for E' z is the root of 
equation 1761 that lies between —qi/pi and —q e /p e . 



8 SHOCK CALCULATIONS 

All fluid variables are made dimensionless by expressing the 
velocity, magnetic field, electric field, and pressure in units 
of v s , Bo, VsBo/c, and Po/{pv1), r espectively. T he charac- 
teristic length scale of the shock is JWardleHl99sh 



VAO 



Equation 1761 h as N poles at E' z = —qj/pj, with N— 1 solu- accuracy is reached. Once E' z and /3i, e , g are known, v;, e , 9 — v 

are calculated (equation 1191 ). enabling the calculation of 
the current density (equation 1131 ). The heating and cool- 
ing rates, (equations 1451 and (15010 . can also be found. The 
right hand sides of equations 1281 . ©, and J7J, can then be 
evaluated and passed to the integrator. 

The calculation of radiative shocks is complicated by 
the long length scales required over which the shock cools 
to its post-shock temperature. The jump conditions yield 
the downstream state of the shock but integration cannot 
proceed directly from this state as the length scales for the 
shock to cool after the other dynamic variables relax is too 
large. Relaxation m ethods have been implemented in the 
past for this reason (|Draini|T980), requiring the determina- 
tion of the fluid variables on grid points in B x — B y — P space 
and then the variables are 'relaxed' until the correct solu- 
tion is reached that satisfies the boundary jump conditions. 
This method is inapplicable here since the calculation of E' z 
requires an initial guess at each point on the shock and an 
arbitrary choice of E' z will not always satisfy bracketing pole 
criteria of equation 1761 . Integration through the shock how- 
ever, provides the initial guess for E' z , which is sufficiently 
close to the desired solution and satisfies the pole criteria. 

The alternative we have implemented is to search for 
the solution from upstream of the final post-shock state, be- 
fore radiative cooling is dominant. An initial starting state 
for (B*, B*, Pd") is chosen in which B x and B y should 
have almost obtained their post-shock values; B y m and 
B* = B xd + 5B X , where 5B X is small ( ~ 1% of B xd ). The 
pressure P* — Pd + SPd is given a larger relative pertur- 
bation and the radiative cooling zone is stepped over. It is 
unlikely that integration from this approximate state will 
yield the correct solution, but it acts as a good initial guess. 
The calculation is a two point boundary problem with a 
known desired final (upstream) state and the (downstream) 
initial values need to be adjusted until integration from 
downstream to upstream yields the correct fast shock solu- 
tion. This is done by using the shooting method (iPress et alJ 
1992) which repeatedly performs the integration for adjusted 
initial states. The initial starting parameters (B x * , B y * , P£) 
are free variables, and the upstream variables (B x o, B y o, Po) 
are fixed. The shooting method proceeds by using a multi- 
dimensional Newton- Raphson root finder which zeros three 
functions obtained by integrating the ODEs in B x , B y , and 
P over the domain of integration. 

A requirement of the shooting method is that for each 
initial guess (B x , B y , Pd), the subsequent shot can trans- 
verse the entire domain of integration. This is not always 
the case for very wrong initial guesses, and physically im- 
possible states may be reached (e.g., T or x e may go neg- 
ative). The convergence of the solution is complicated by 
unsuccessful integration shots. Fortunately, for any number 
integration shots there was clear convergence, before the tra- 
jectories diverged. A multiple shooting method was therefore 
implemented which after convergence did occur, restarts the 
shooting process at a point inside the shock using values de- 
rived from the converged solution from the previous shooting 
procedure. This process is repeated through the domain of 
integration, until the upstream state is reached. 



E,-7joPjo' 



(79) 



The shock models are calculated as follows. For the MRN 
grain size distribution models, x g o is assigned to each grain 
class (Section 15.21 . x e o is then calculated from charge neu- 
trality (equation 1141 ). and /3 e o, /3io and (3 gm o are calcu- 
lated (equation I2H ). The appropriate root of the isother- 
mal jump condition (equation 1401 ) for downstream v z d is 
found, allowing the calculation of downstream v x P, and 
B x from equations 1251 . 1271 and 139H . The pole bracketing 
species for E' z is ide ntified a s desc ribed in Section |7| 

As discussed in IWardld (1998), a one parameter family 
of solutions exist and integration from upstream to down- 
stream yields the intermediate shock solutions, and from 
downstream to upstream gives the desired fast shock solu- 
tion. Integration of the ODE's for B x , B y , and P, equa- 
tions ©,0, and 1281 1 . respectively, are then performed af- 
ter perturbing these quantities from the downstream state. 
The e quations are stiff so the integration method of iGearl 
lll97lf) is used. For each step, the integrator produces val- 
ues for B x , B y and P, enabling the calculation of v x , v y , 
v z , E' x and E' y directly from equations I25t - I27j and (1321 - 
13311 . Finding solutions to equation 1701 for f} e and /3 g k and 
equation 176K for E' z , are complicated by their mutual de- 
pendence. To find a solution, the value of E' z obtained in 
the last step of the integrator is used as a starting value. 
The corresponding Hall parameter (equation 1691 ) is found 
(using the updated values of B x , B y , E' x , and E'y), including 
the corresponding bracketing poles for equation 1761 . The 
accuracy of this value for E' z as a root of equation 176K can 
then be tested. If tolerances are not met, the value of E' z is 
amended and the process repeated until required numerical 
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z (10' 6 cm) 

Figure 1. Single sized grain shock model with grain radius 0.1/xm, n g /riu = 1.6 X 10~ 12 , and 
components B x , B y , the neutral pressure P, ratio c s /v z and temperatures T,Ti and T e . Also plotted are the electric field components 
in the neutral frame E' x , E' y , and E' z , the neutral velocity components v x ,y,z, and the charged species drift velocity components 
(vix, y, z — Vx,y,z) and (v gXy y iZ — V x ,y,z)- The grain parameters Z g and f3 g , are shown along with the compression ratios v B /vj z and the 
heating rates per unit mass, G nJ /p. 
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9 RESULTS: C-TYPE SHOCK MODELS 

C-type shock profiles are presented assuming both small and 
large single sized grains (Sections 19,11 and 19. 3^ as well as 
the MRN grain size distribution with and without man- 
tles and PAHs (Sections E31 is varied in Sec- 
tions 19.21 and 19.41 and the effects of suppressing B y are 
demonstrated in Section 19.71 Unless stated otherwise, the 
pre-shock conditions are: Bo = 0.3mG, nu = 10 5 cm~ 3 , 
n = n(H 2 ) + n(H) + n(H e ) « 0.5n H + 0.1n H = 0.6n H , 
nm/nH — 3 x 10 -8 (by equation 1521 1. 7 = 7/5 (see the 
discussion near equation llll l. and To = T e o = To = 10K. 
The shock speed is v s — 18km/s 



9.1 Single sized small a g = 0.1/im grain model 

A single grain species, a g — 0.1/i,m and internal den- 
sity of 2.5g/cm 3 , is now considered. If the total mass o f 
the grains is 0.01 that of hydrogen jDraine et alj fl9831. 
n g0 /n H = 1.6 x 10" 12 JWardlelll998Tl . n i0 /n H = 3 x 10" 8 
and n e o/riH = 2.99984 x 10 -8 by charge neutrality (equa- 
tion dj). Then (3 i0 = 1.382 x 10 4 , f3 e0 = -3.358 x 10 7 , and 
/3 9 o = -1.632. Fig. □ shows a shock profile with 9 = 40°. 
The B x — By phase plot shows the shock transition from the 
upstream state (B x o,B y o) = (0.7, 0)-Bo to the downstream 
state (B x d, B y d) = (13.7, 0)-Bo- The rotation of Bx out of 
and back into the x — z plane is easily identified. 

Initially B y then B x increase (from upstream to down- 
stream) near z — 1 x 10 16 cm, resulting in an increase of B±. 
As B\/(2Bq M\ ) +v z /v s + P/ (poVs) is constant (by equa- 
tion I27H . and changes in P/(pov 3 ) are relatively small in 
this part of the shock there is a compensatory decrease in 
v z /v 3 . The compression of the neutral fluid begins, and P 
and T begin to increase. T e and T also increase, and the ap- 
proximation used for T e (equa tion I|55|l1 leads to a maximum 
of ~ 4000 K, consistent with lDraine et alJ dl983ll . 

The ODE for P is sensitive to the compression of the 
neutral fluid, and the balance of G — A with the magnetic 
gradient terms of equation 1281 . The G — A term dominates 
in equation 1281 . Initially P increases as the neutral fluid is 
compressed and heated. T decreases after z = 2.05 x 10 16 cm, 
but P continues to rise since there is still compression. P 
does not fall until z — 2.4 x 10 16 cm when the compression 
slows, and v z /v 3 reaches its downstream value. T and P 
have not yet reached their post-shock states in Fig. and 
will not reach them until far downstream of the shock, after 
the neutral fluid cools inside the radiative zone. 

The P profile shows two discontinuities in dP/dz at 
z = 2 x 10 16 cm and z = 2.1 X 10 16 cm, indicated by dotted 
lines. These occur at T — 1087K and are artefacts of the 
cooling functi on (equatio n 1)501 1 which is discontinuous at 
T = 1087K jLepp fc Shulll ll983l. Interpolation of A was 
used to smooth the transition at T = 1087K, however P is 
highly sensitive to A, still reflecting the discontinuity in A. 

Fig-0shows the ratio c a /v z nearing the sonic condition 
at z — 2.45 x 10 16 cm (c s /v z ~ 0.9). The peak coincides with 
the slowing of the compression as v z approaches the down- 
stream value. c s /v z then decreases as c 3 decreases with T. 
v z = c s if the neutral fluid cools quickly enough after the 
peak in T. This occurs for certain pre-shock conditions, e.g., 
high pre-shock number densities and/or low 0. In this case, 
the shock becomes either C*-type or J-type, requiring dif- 



ferent numerical procedures to calculate the shock structure 
since equation 1281 c annot be int egrated through the critical 
point jRoberee fc Draind Jl990M . 

In FigQ E' x and E' z are similar in magnitude, with E' y 
an order of magnitude larger. Since v x oc (B x — Bq x ) and 
v y tx By (equation 1261 1. trends in v x and v y follow those of 
B x and B y . v x steadily increases through the shock, before 
reaching the downstream state given by the jump conditions. 
The magnitude of v y is orders of magnitude smaller than v x 
and v z . As E'y = (B • E')B/P 3 is negligible through the 
shock, and by equation 1191 



0J (E'xB) 



B- 



+ 



1 + 0] B ' 



(80) 



For \/3j\ » 1 the species are tied to B with (vj — v) w 
(E' x B)c/B 2 , as is the case for the ions and electrons. The 
ions (and electrons) lead the neutrals in the x and y direc- 
tions but lag the neutrals in z. The second term in equation 
1801 becomes non- negligible for the grains with \/3 g \ ~ 1 — 10. 



(v iy 



has more structure than both (yi X 



The 



shoulder at z ~ 2.2 x 10 16 cm where (vi y — v y ) declines much 
less rapidly, corresponds to the peak of E' z and a slowing 
down of the compression of B x (B x ~ B x d)- As E' z subse- 
quently decreases after z = 2.25 x 10 16 cm, there is a corre- 
sponding reduction in (vi y — v y ). 

are shown in Fig. and are similar in 



magnitude to (vi: 



is an order of mag- 



nitude smaller then (vi y — v y ) and contains oscillations. 
The ratios v s /vj Z (equations 1221 and 1231 1 show that the 
grains and ions are compressed ahead of the neutrals at 
z = 1.8 x 10 16 cm (along with B), with compression in the 
neutral fluid occurring after z — 2.0 x 10 16 cm. 

The grain charge remains —1 until T e increases and 
\Z g \ > 1 by equation 1561 . The peak in \Z g \ corresponds to 
the peak in T e . Initially \(3 g \ ~ 1 , but as the grain drift in- 
creases, u g increases (equation 1171 1 and there is a decrease 
in \P g \. \P g \ starts to increase when \Z g \ increases. As B is 
compressed \(3 g \ increases further and the grains become bet- 
ter coupled to B. The decrease in \(3 g \ after z = 2.1 x 10 16 cm 



follows the decrease in \Z„ 



AS \{ig\ 



10, v 



as T e declines. 

- v is dependent on the balance 



of terms in equation 1801 . When |/3 S | w 1 the second term 
reaches a maximum for a given B and E', and v 9 — v is 
determined by the balance of both terms. For |/3 9 | < 1 both 
terms decrease with \j3 g \ and |v 8 — v| decreases, signifying 
an increase in the grain-neutral coupling. The oscillations in 
(v gy — v y ) axe due to the balance of the competitive terms in 
equation 1801 . The sharp increase in (v gy — v y ), marked by 
the dotted line, corresponds to the increase in \Z g \ and |/3 9 |. 
G" s dominates G in Fig. with only a small contribution 
from G ni . 



9.2 a g = 0.1/im grain model: varying 6 

A comparison of the a g = O.lum models for 8 — 
40°, 45°, 60°, and 9 = 80°, is plotted in Fig.0 The 2 offset 
in each 6 model is arbitrary, and the solutions have been 
shifted in z to overlay and compare. Since B x o = Bq sin 9 
the upstream state differs for each model. As 9 decreases 
from 80° to 40°, the amount of Bx rotation increases as 
seen in the B x — B y phase plot.. For the 9 = 80° model, B y 
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Figure 2. Comparison of 9 = 40° (solid), 45° (dot-dash), 60° (dotted), and 80° (dashed) for the a g = 0.1/xm and n g /n H = 1-6 X 10" 
model. The 8 = 40° line series is the same as in Fig.Q The descriptions of the quantities are the same as in the caption for Fig. Ill 



and v v are small, and the neutrals are restricted mainly to 
the x — z plane. Conversely, for 6 = 40° there is maximum 
rotation of Bi with largest B y and v y . The shock widths do 
not differ greatly with 6. There are larger peaks in P and 
T as 8 is decreased. The neutrals approach a sonic point as 
9 decreases, with the peak c s /v z — * 1. C-type solutions are 
restricted by the choice of 9, 9 crit i Ca i ~ 40°, below which the 
neutral fluid becomes subsonic inside the shock. 

The jump conditions are 9 dependent so the down- 
stream values of v x and v z differ for each 9 model in Fig. |2] 
The differences in B y and B z with 9 leads to differences in 
(vix ~ Vx), e.g., when 9 = 80°, the peaks in B z and B y , and 
thus (E' x B)^ are small, leading to a lower peak (vi x — v x ) 
(equation 1801 1 ). (vi y ,z —Vy lZ ) show smaller differences, with 
the largest peak in \v%y, z — v y,A for the low 9 — 40° case. 

(v g x,z — Wi, z ) are similar in magnitude to (vi XiZ — v XiZ ), 
showing the same 9 dependent features. Differences in the 
structure of (v gy — v v ) with 9 are large, and as 9 increases 



from 40° to 80° the number of oscillations decreases and the 
relative peak magnitude increases. \Z g \ and \/3 g \ axe compa- 
rable for all 9 (see Fig.0, and differences in (v gXlV]Z — v X]ytZ ) 
are mainly due to the differences in E' and B. 

9.3 Single sized large a g = 0.4/^m grain model 

Dust grains inside dense molecular clouds are larger than 
those in the diffuse interstellar medium (e.g.. lCarrasco et alJ 
Il973l) due to coagulation and the growth of ice mantles. 
The grain size is now increased to a g = 0.4/xm to represent 
grain growth via coagulation. Assuming the total mass of 
the grains is 0.01 th at of hydrogen then x g o — 2.5 x 10 -14 
(Drai ne et al.lll983i) and /3 9 o = —0.1020. x g o and the net 
grain charge is decreased by two orders of magnitude com- 
pared with the a g — 0.1/im model. The a g = O.lfim and 
a g — 0.4/^m models for 9 = 45° are compared in Fig. [3] The 
0.4/im grains are less coupled to B, so the rotation of Bi 
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increases. The shock width is larger and there is a corre- 
sponding reduction in the peak temperature to T ~ 800K, 
as the same amount of kinetic energy flux must be dissipated 
into thermal energy as for the thinner 0.1/im model shock. 

The remaining variables plotted in Fig.E]are consistent 
with the field components, e.g., there are larger peaks in 
\ v iy,gy ~ v v\ an d v y, consistent with larger B y . There is also 
a steeper decrease (increases) in v z (v x ) for the smaller grain 
case then for the a g = 0.4/im case. The maximum peaks of 
(vix,z — v x , z ) and (v gx , z — v x ,z) are similar for both models. 

The decoupling of the 0.4/im grains is clearly seen in the 
ratio v s /v gz (equation 12311 ) . Initially the grains are com- 



pressed along with B and the ion fluid {v s /v g 



e/Viz) 



until after z~3x 10 16 cm when |/3 9 | and \Z g \ decrease. In 
which case, the grains are decoupled from B and the com- 
pression slows compared with the ions. The grain fluid un- 
dergoes expansion until z = 3.7 x 10 16 cm. Compression then 
begins again, shown by increasing v s /v gz , with the transi- 




10.0 



o.(10- 5 cm) 



Figure 4. Grain radius a g plotted with x g o = n g o/nu for a stan- 
dard MRN distribution (equation I6UI ) and a MRN distribution 
with grain coagulation and ice mantles using 10 grain size classes. 
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Shock profile with 9 = 45° and the MRN grain size distribution of Fig. [I] The descriptions of the quantities are the same as 
The grain parameters for each size class are plotted in a different colour, with the radius a g indicated in the (v gz — v z ) plot. 



tion of \/3 g \ from |/3 9 | > 1 to \/3 g \ < 1 as the grains become 
completely decoupled from B, and compressed along with 
the neutrals after z = 3.7 x 10 16 cm. 

A lower grain abundance for the 0.4^im model leads to 
a reduced heating rate per unit mass G" 9 , whereas G nl are 
similar for both models. The neutral temperature T (and 
Ti >e ), reaches smaller peak values, however the 0.4/im grains 
obtain larger charge \Z g (a, g ,T e )\ as they have a larger surface 
area for electrons to stick to (equation 

9.4 MRN grain size distribution 

In this section a MRN grain size distribution is adopted 
with 10 grain size classes. x g o — n g o/riH versus size class a g , 
(equation IbUH 'l is plotted in Fig. 3] for the standard MRN 
grain size distribution. The total pre-shock grain density is 
J2 g n s0 = 3.39 x 10" 5 cm" 3 , leading to x e0 = x i0 -J2 g %g0 = 
2.966 x 10 -8 (equation p4t ). In general, the smaller a g , the 
larger x g o- The only exception is the 65 A class with x g o 



greater than that of the 53A class, since for smaller a g the 
Gauss-Legendre weight w m is less than the corresponding 
W m for midrange a g (equation I6UH ). A significant propor- 
tional (67%) of the total grain number density n g o comes 
from the smaller a g S~ 65A grains, which couple best to the 
magnetic field and suppress B diffusion. As the number of 
grain bins is increased from 1 to 10, the shock solution con- 
verges; The shock solution is already converged for > 10 
grain bins (plots not shown here). 

Fig- E]shows the corresponding shock profile. The max- 
imum peak in B y is significantly smaller than that of the 
equivalent 0.1/im grain model (Fig.|5J. The larger (low |/3 9 |) 
size classes are better coupled to the neutrals through col- 
lisions, but their abundances are too low to have any large 
effect on the overall structure. The small dominant grain size 
classes (large \f3 g \) remain coupled to B and so there is min- 
imal rotation of B out of the x — z plane, and subsequently 
a lower peak in B y . The shock is also thinner and the peak 
in T is higher, compared with the single grain models. 



Dust grain dynamics in C-Type shock waves in molecular clouds 13 



1500 







1^ 






/ 1 






/ 1 


\ j 
\ \1 









0.00 



0.20 0.30 0.40 
z (10 ,6 cm) 



0.00 




0.20 0.30 0.40 
z (10 16 cm) 



0.20 0.30 0.40 
z (10 ,6 cm) 



i.50 0.60 



Figure 6. Shock profiles for the MRN grain size distribution (solid line) and MRN a distribution with mantles (dashed line), for 6 = 60° 

same as 



and x g o from Fig.|H (vgz,y,z 
in the caption for Fig. [F] 



■ t z) are shown for the MRN(mantles) model only. The descriptions of the quantities are the ; 



In Fig. |S] (v gx — t) x ) clearly differs with a g , and the 
smallest grains are accelerated ahead of the larger grains 
until the peak at z = 0.63 x 10 16 cm. The \(3 g \ plot shows 
for the smallest grains \/3 g \ » 1 through the entire shock 
with (v gx — v x ) ~ (vix — v x ). Conversely, for mid size range 
grains with \/3 g \ « 1 initially, |/3 9 | reaches a maximum of 
only ~ 10 due to increases in \Z g \. The largest grains have 
the greatest surface areas, and so obtain largest \Z g \. As 
\/3 g \ increases, the larger grains become better coupled to 

B and (v gx Vx^larae grain * i^gx ^x) small grain- At Z = 

0.63 x 10 16 cm, (v gx — v x ) is a maximum for all sizes, at which 
point v x increases, and (v gx — Vx) subsequently decreases, 
and all grains classes are (near) uniformly decelerated in x 
with respect to the neutrals. Similar comments can be made 
in regard to the magnitudes of (v gz —v z ). 

(vgy — Vy) AiRer subst ant ially with a g in both sign and 
magnitude. The smallest a g {\j3 g \ » 1) class have (v gy — 
v v) ~ ( v iy~ v y) as expected, and are accelerated ahead of the 
neutrals. However the largest a g (\/3 g \ ~ 1 — 10) class obtains 
large negative (v gy —v y ) (as was also observed for the single 

0. 4/im model). For each of the larger grains a g > 825A there 
is a sharp trough in (v gy — v y ) which coincides with \/3 g \ ~ 

1, after which as \/3 g \ increases, the grains become better 
coupled to B and \v gy — v y \ decreases. 

Each grain class is differentially compressed as shown in 
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the v s /v gz plot (equation 1231 ) in Fig. [S] The compression 
of the larger a g classes begins to slow (d(v z ~ 1 )/dz ~ 0) near 
z = 0.8 x 10 16 cm, before increasing as v gz — » v za - - Initially as 
\P g \ increases, the larger grain classes become better coupled 
to B and are compressed along with the ions and B. The 
compression then slows with decreasing \(5 g \ after z = 0.7 x 
10 16 cm, and the larger grains decouple from B. As f3 g ~ 1 
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the compression slows (d(v~ 1 )/dz ~ at z = 0.8 x 10 16 cm), 
after which |/3 9 | decreases further and the large grain classes 
are compressed along with the neutrals (see also Fig. 19.31 . 
In general, the peak in G Ng / p, also shown in Fig.|S] increase 
with decreasing a g (due to increasing x g ), and the smallest 
grains dominate the heating. G nl is comparable to that of 
the a g — 1920A class. 

9.5 MRN grain size distribution with mantles 

The MRN distribution with ice mantles (discussed near 
equation 1591 ) is now considered. x g o for each size class 
are shown in Fig. 0] with J2 g n g0 = 3.39 x 10~ 5 cm" 3 . A 
comparison of shock profiles for 9 — 60°, using MRN distri- 
butions, with and without grain mantles are shown in Fig. 
|S| The rotation of Bi is increased for the MRN(mantles) 
case as the grains size classes are larger and are only par- 
tially coupled to B. The shock width is narrower with a 
corresponding higher peak T. A higher peak in T leads to 
higher c 3 , and shock solutions are limited to 8 cr it ~ 60° for 
the MRN (mantles) model, as opposed to 8 cr it ~ 45° for the 
MRN case. A larger peak in B y leads to higher peaks in v y 
and (v iy — v y ). (v ix , z — v x ,z) and v x ,z are similar in mag- 
nitude for both models. In general the grain drift velocities 
show similar trends to the MRN model (Fig. 01. 



The B x - By phase plots for the MRN and 
MRN (mantles) models with varying 9 are compared in Fig. 
[7| In general, there is a decrease in the rotation of Bx with 
increasing 9. The corresponding effects on the neutral and 
charged species dynamics follow on from the relative mag- 
nitude of By through the shock. MRN (mantles) 9 — 60° is 
essentially equivalent to MRN 9 — 45° . Increasing the grain 
sizes, which leads to larger forces out of the shock plane 
and a more non-coplanar solution, is essentially equivalent 
to taking a MRN distribution with smaller 9. This result is 
initially surprising as one might expect the addition of more 
highly coupled charge carriers with \/3 g \ » 1 to reduce the 
Hall current. However electron abundance and therefore net 
charging of grains is reduced, so many grains remain decou- 
pled from B. 



9.6 MRN with PAHs 

In the interstellar medium, for ng = 10 5 cm -3 , the abun- 
dance of negatively charged PAHs" is of the order 10~ 8 , 
with a similar total abundance, 1.1 x 10~ 8 , of ions and 
positively charged PAH + s (values taken from Fig. 1 of 
iKaufman fc Neufeldl Jl99rj) ). Assuming a MRN distribution 
as before, ^2 g x g o = 3.39 x 10 -10 , a PAH - component; 
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x V aho — 1 x 10~ 8 , and the corresponding PAH + and ion com- 
ponent; Xio = 1.1 x 10 -8 , then by cha rge neutrality, x e p = 
6.606 x 1CP 10 (consistent with x e o in iKaufman fc Neufeldl 
( 1996)). Since x e o is two orders of magnitude lower than in 
the previous models ( Sections 19 . 119 . 5|l . there are not enough 
electrons to charge all grain classes inside the shock. Even 
if only the smaller grain size classes are charged, which ac- 
count for a significant fraction of the total grain surface area, 
these are already highly coupled to B (\/3 g \ » 1) and in- 
creasing \Z g \ does not effect the grain dynamics or the shock 
structure. For simplicity, Z g is assumed constant (Z g = — 1). 
The rate coefficient for PAH-neutral scattering is t he same 
as for the ions with < av > pah = 1.6 x 10 cm s" - 1 iWardld 
I199ST) . and /3 pah0 = -fto = -1-382 x 10 4 . 

The shock profiles (6 = 45°) for MRN(PAH) and MRN 
are compared in Fig.|H| There is significantly more Bi rota- 
tion for the MRN (PAH) model. The B x and T profiles are 
similar, leading to the close similarities in v x (v% x —v x ) and 
v z (viz —Vz). Ions, electrons, and PAHs, remain well coupled 
to B , and (v< - v) = (v e - v) = (v pQh - v). 

There is a clear distinction between {v gXtZ — v x ,z) f° r 
the smaller and larger grain classes. The smaller grains re- 
main well coupled to B, whereas, the larger grains remain 
highly coupled to the neutrals through collisions and are co- 
moving with them (v gx , z — v x ,z ~ 0). With fixed Z g , the 
larger grains do not become increasingly coupled to B as 



they pass through the shock (as in the previous results for 
MRN). From upstream to downstream, each smaller grain 
is compressed ahead of the next larger grain. The small- 
est grains are compressed with the ions and magnetic field 
through the entire shock, and the larger grains are com- 
pressed with the neutrals. The midrange size class ratios 
Vs/vzg he between these two extremes. 

(v gy — v y ) show (in general) similar relationships be- 
tween the relative magnitude and sign with grain class size, 
as seen previously in Fig.0 Smaller grain classes are coupled 
to B with (v gy — v y ) ~ (viy — Vy) > 0, and the larger grain 
classes have (v gy — v v ) < 0. There is also a turnover in the 
magnitude of (v gy —v y ) with increasing a g . The larger grain 
size classes (|/3 S | << 1) are highly coupled to the neutrals 
and have small \v gy — v y \. As a g decreases (\/3 g \ ~ 1 — 10), Hall 
effects are present and these obtain larger negative {v gy —v y ). 
As a g is decreased further, the coupling to B increases, and 

{Vgy - Vy) -» (v iy - Vy). 



9.7 Suppressing B y 

Fig.|51shows the effect of suppressing B y in the MRN (PAH) 
model. v x ,y,x, (v gx ,z —v x , z ) and (v ix ,z —v XlZ ) are essentially 
unaltered, except for the slightly higher peak in (vi X ~ v x ) for 
By — 0, due to the corresponding differences in (E' x H) x . 
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The peak in T is lower for B y — 0, but only by ~ 10K. 
The heating rate is increased, however the cooling rate is 
a strong function of T and so the change in T is small. 
In y all grain size classes lag the neutrals (v gy — v y ) < 0, 
and for the smallest grains (v gy — v y ) ~ 0. The larger grain 
size classes have similar (v gy — v y ) as in Fig. [8] as they are 
decoupled from B, and so are unaffected by changes in B y . 
The corresponding B y — results of the models in Sections 
l9.H9.6l plots not shown) also show only minor differences. 



10 DISCUSSION 

The results above demonstrate the effects of charged dust 
grains on the structure of fast C-type shocks. For uh = 
10 5 cm~ 3 and shock speed v s = 18km/s radiative C-type 
shock solutions exist {8 critical < 8 < 90°). The grain model 
details, i.e., modelling the grain population as either a sep- 
arate grain species or else as a continuous grain size distri- 
bution, is important to the overall shock structure, as well 
as the dynamics of the individual grain species/size classes. 

The grain models considered here were chosen to illus- 
trate the main qualitative changes in the resulting shocks. 
The shock profiles were clearly altered with the introduction 
of a grain size distribution as opposed to the single sized 
grain models. The net charge on the grains is increased for 
the MRN models with the grains becoming better coupled 
to the magnetic field, resulting in less magnetic field rotation 
inside the shock front. The MRN (mantles) and MRN (PAH) 
models were chosen to illustrate changes in the shock profiles 
when departing from the standard MRN model. It should be 
noted that the MRN (mantles) model is not purely a MRN 
model with growth of ice mantles of constant thickness. 
We have followed the model o flNishi. Nakano fc Umebavashil 
il99lf) who increase the size of the grains by a factor of 9/5. 
The qualitative changes are clear though; any increase in the 
grain sizes leading to larger grain surface areas (larger grain 
charge) will cause an increase in the net charge on the grains, 
and the grain coupling to B. However, by increasing the 
collisional cross section, the grains are increasingly coupled 
to the neutrals via collisions. In the MRN (mantles) models 
shown in section l9~51 the latter effect was dominant with the 
MRN (mantles) solution showing more magnetic field rota- 
tion than it MRN counterpart. 

As B is compressed inside the shock, the charged species 
are accelerated with respect to the neutrals, however, this is 
opposed by the collisional drag of the neutrals. If the grains 
are highly coupled to B (high \f3 g \) and the movements of the 
neutrals and other highly coupled species (i.e., ions, PAHs, 
and electrons) are mainly restricted to the x — z plane, then 
B y is small. If however, the grains are only partially coupled 
to B, either because they are large in size and so are better 
coupled to the neutral fluid via collisions and/or have low 
charge \Z g \ (and low \/3 g \), then the grains drift out of the x — 
z plane, and the shock becomes non-coplanar with large B y . 
The grains dominate the collisional heating of the neutrals, 
with only a small contribution from the ions and a negligible 
component from the electrons (due to their low mass). 

The non-coplanarity in the shock solution depends on 
the grain model; the single grain models (Sections I9.ll9.3t 
produce significantly more rotation of Bi than for the 
MRN models ( Sections 19". 419.61 . For both 0.1/im and 0.4/um 



models, the grains were only partially coupled to B with 
\/3 g \ ~ 0.1 — 10 inside the shock and the out-of-plane forces 
are substantial. For the MRN models, rotation of B is mini- 
mal, since the more abundant, smaller a g classes {\(3 g \ » 1) 
which are highly coupled to B, dominate the net grain ef- 
fects on the shock structure. 

As 8 — > 90° the peak in B y decreases, for a given grain 
model, and the well B coupled charged species are restricted 
to the shock plane (weakly coupled grains are still able to 
drift out of the shock plane). Bi rotation increases with de- 
creasing 8 and |vj — v become larger (on average), resulting 
in more dissipation. The shock must dissipate the incoming 
kinetic energy flux over a shorter length scale, and the peak 
temperature T inside low 8 shocks (for a given grain size 
model) are therefore higher than for high 8 shocks. For low 
8 shocks, the neutral fluid experiences more compression as 
it cools and so the sound speed is higher in the final part 
of the flow. There is a critical point for low enough 8, that 
in the neutral fluid v z — c s , in which case the derivative 
dP/dz is ill-defined (equation 1281 ') and the shock solution 
is no longer C-type. Below this critical angle 8 < 8 cr iticai 
the shock solutions may be either J-type or C*-type. The 
limiting value of 8 cr iticai varied with grain model, but the 
ordering remains same, the higher the peak temperature in 
the shock, the larger 8 critical- 

Grains dominate the collisional heating of the neu- 
trals. The peak T was higher and shock width narrower for 
the 0.1pm model (Section 19. II compared with the 0.4pm 
model (Section 19.31 since the pre-shock abundance of the 
0.1/im grains was two orders of magnitude larger, than for 
the 0.4/im model, which compensated for their smaller col- 
lisional cross-sections. Conversely, for the MRN (mantles) 
model fSection !9.5l . the increase in the total frictional heat- 
ing of the neutrals due to the increase in the net grain col- 
lisional cross-sections lead to hotter thinner shocks than for 
the MRN model with an equivalent pre-shock grain abun- 
dance (Section 19.41 . Since the temperatures, and subse- 
quently c s , are higher in the MRN (mantles) model (for a 
given 8), it is more difficult for the neutral fluid to main- 
tain supersonic flow and subsequently 8 is restricted to 
8 > 8 critical ~ 60° for the MRN (mantles) model but only 
8 > dcriticai ~ 45° for the MRN model (for v 3 and uh con- 
sidered here). 

There are also critical values of uh for which C-type so- 
lutions do not exist. Increasing nn (for a given 8), increases 
the neutral-grain collisional heating rates and thus dissipa- 
tion. With a large enough increase in T and c s the neutral 
flow cannot remain supersonic. For all grain models consid- 
ered here, the C-type solutions also broke down for uh > 
10 5 cm -3 (independent of 8), with only J-type or C*type 
shocks possible. The exact limits of 8 and uh could vary 
with the choice of cooling function. The lLepp fc Shulll ll983) 
cooling rate has been used here, other treatments of the 
coolin g rate dMartin. Schwart fc Mandvll996l ; lGalli fc Pallal 
1998), will lead to different temperature structures and thus 
different limits on 8 and uh- However the orderings of the 
different grain models will remain the same; the cooler MRN 
and 0.4pm models will have lower 8 cr iticai than the hotter 
MRN(mantles) and 0.1pm models, respectively. 

In the MRN models (Sections 19.419.51 the T e depen- 
dent grain charging, Z g (T e ), increased the grain coupling to 
B. The smaller a g classes \(5 g \ » 1 were compressed with 
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the ion fluid and B. Mid range grains remained coupled to 
B until reductions in Z g lead to their decoupling from B, 
and were compressed with the neutral fluid thereafter. The 
largest grains \/3 g \ « 1 are essentially compressed along 
with the neutral fluid. When PAHs are present fSection l9.6H . 
x e o is reduced with Z g = — 1, and the mid to large grain size 
classes are less coupled to B and are (on average) decoupled 
further by neutral collisions inside the shock. As a result 
there is a larger Hall current and the amount of rotation of 
Bi is increased. 

Here T e is given b y the approximation in equation 15511 
based upon Fig. 1 of iDraine et al.l (11983ft . Changing the 
treatment of T e could largely effect Z g of the larger grains 
and alter their coupling to B as well as (v 9 — v). However, 
smaller grains with low Z g are highly coupled to B (essen- 
tially independent of T e and Z g ) and dominate the grain 
effects on shock structure. Thus variations in the treatment 
of T e are unlikely to change the model results. 

\Z g \ increases with increasing T e according to equation 
15611 . however as T e decreases, the uncharging proceeds sim- 
ply via the re-release of the captured electrons. More realis- 
tically the reduction in \Z g \ proceeds via the capture of ions 
iDraine fc Sutin|ll987f) . and the ion fluid is depleted. More- 
over, inside shock waves in dense clouds, negatively charged 
PAHs may be neutralized in neutral-PAH~ reactions, result - 
ing in electron detachment iPineau des Forets et alJll988T) . 
iFlower fc Pineau des Forets! ]2003ft" modelled the charge 
variations of grains in a planar steady C-type shock with 
v 3 = 50km/s. The grain population was represented by a 
small (dominant) PAH grain component and a large grain 
component (given by an MRN distribution 0.01 ^ a g ^ 
0.3/im), and they found that the charge variations for each 
were very different. The PAHs were neutralised by neutral- 
PAH~ reactions via electron detachment, while the larger 
grains became more negatively charged as T e and thus \Z g \ 
increased through the shock. 

Detachment of electrons via neutralization of the 
PAH~s, may allow for the Z g (a, T e ) (equation 15611 ) charging 
to proceed in the MRN (PAH) model, potentially increasing 
the grain-B coupling, reducing the rotation of B^. However, 
the capture of ions cou ld also counteract any increases i n 
Z g . The shock model of lFlower fc Pineau des Forets! J2003D . 
included both electron and ion attachment on the grain sur- 
faces, and the net negative charge on the grain component 
still increased significantly due to the increase in the attach- 
ment rate of electrons with increasing T e . Thus, it is likely 
that if ion attachment were included in the models here, \Z g \ 
should still increase with T e . The details of this will need to 
be quantified and confirmed in further studies. 

Inertia of the grains has also been neglected and their 
drift is given by the balance of electromagnetic forces 
and neutral-grain collisional drag. lCiolek fc Robergel (12002ft 
demonstrated that for very small grains (a g < 0.1/im) the 
length scales over which the neutral fluid and magnetic field 
were compressed, L neut and Lf ie i d , respectively, are over 
three orders of magnitude larger than the length scales for 
the gas drag to decelerate the grains Ld rag (g)- Thus, neglect- 
ing the inertia of the smaller grains is a good approximation. 
However, for the large grains (a g ~ 0.1/^m), Ld rag {g) is of the 
same order as L f ie id and so the inertia cannot be neglected. 
The larger grains dominate the field-neutral coupling and so 
the inclusion of their inertia may effect the shock structure. 



The amount of dissipation inside the shock was rela- 
tively insensitive to the orientation of Bo for the pre-shock 
conditions considered here; tih = 10 5 cm -3 . The peak T 
inside the shock was essentially unaltered when B y = 
for the MRN (PAH) model (Section |53 even though the 
MRN(PAH) model had the most Bi rotation. Similar com- 
parisons have been made for the other MRN models. The 
critical angle 8 cr iticai for which the neutral flow inside the 
shock makes the transition from supersonic to subsonic is 
also unaltered when B y is suppressed. For the purposes 
of calculating detailed chemical models and molecular line 
emission the shock problem can be greatly simplified by sup- 
pressing By. This may not be the case, however, for higher 
uh where grains may become further decoupled from B. 

The choice of grain model affects the shock structure, 
however, the same overall jump conditions (for given pre- 
shock conditions) are still satisfied, and the same power per 
unit area is dissipated inside the layer of hot gas inside the 
shock transition. For the models here, the temperatures in- 
side the shock reach ~ 1500K and the molecules will remain 
intact, so the molecular lin e emission should be sim ilar re- 
gardless of the grain model. iNeufeld fc Stond (119971) numer- 
ically modelled CO and H2 line emission for C-type shocks 
subject to the Wardle instability, concluding that the line 
strengths are relatively insensitive to the exact details of 
the shock structure. However, the choice of grain model does 
effect the limits ng and 6 for which the shock becomes J- 
type (or C* type), and the temperatures inside J-type shock 
reach much higher values, so the shocked molecules may not 
always remain intact. 



11 SUMMARY 

New shock models for fast oblique C-type shocks have been 
presented here with pre-shock conditions uh = 10 5 cm -3 , 
v a = 18km/s, and Ma = 10 for freely specified magnetic 
field Bo orientation 35° ^ 9 ^ 80°. The grain population is 
represented by either 

• a single grain population of size 0.1/^m or 0.4/im, 

• a continuous (standard) MRN grain size distribution, 

• a MRN distribution with ice mantles, or 

• a MRN distribution with additional populations of pos- 
itively and negatively charged PAHs. 

There were clear differences in the shock profiles with grain 
model due to the coupling of the grains to either/both the 
magnetic field and/or neutral fluid. Smaller and/or highly 
charged grains remain well coupled to B whereas larger 
grains are coupled to the neutral fluid through collisions. 
Shock models with grain populations that are well coupled 
to B result in more co-planar shock profiles with minimal 
rotation of B out of the shock plane. The results are sum- 
marized below: 

• For a given 9, the smaller 0.1 fim grain model produces 
more co-planar solutions than the 0.4/^m model since the 
0.1/xm grains are better coupled to B. 

• Each grain size class in the MRN models were differ- 
entially accelerated inside the shock front, with the largest 
grains being compressed along with the neutrals and the 
smaller grains compressed along with the ions and B. 
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• The more abundant, smaller grains in the MRN distri- 
bution models dominate the net grain effects on the shock 
structure resulting in more co-planar solutions compared 
with the single grain models. 

• The effect of charging grains via the capture of the elec- 
trons increased the grain-B coupling in all grain models. 

• The increase in the grain surface areas from the MRN to 
the MRN(mantles) models increased the net grain charge, 
however the increase in the grain collisional cross sections 
resulted in better coupling of the grains to the neutrals and 
consequently more rotation of B out of the shock plane, for 
a given 6. 

• The grains dominate the drag force and collisional heat- 
ing of the neutrals, so the models with larger (larger colli- 
sional cross sections) or more abundant grains have higher 
peak temperatures and the shock fronts are subsequently 
narrower i.e., the collisional heating for the less abundant 
0.4/im grains lead to lower peak shock temperatures and 
wider shock fronts than for the more abundant 0.1/im grains. 
Conversely, the larger MRN(mantles) model lead to in- 
creased peak temperatures, and narrower shock fronts, com- 
pared with the MRN model, due to the increase in the col- 
lisional cross section of the grains. 

• The higher the peak temperatures reached inside the 
shock front, the higher the sound speed and the neutrals 
cannot always remain supersonic. As a result, the lower 
limit of 6 for which C-type solutions are possible is depen- 
dent on grain model; 9 critical ~ 40° for the 0.1/um model, 
6 'critical ~ 45° for the MRN and MRN(PAHs) models, 
critical ~ 60° for the MRN (mantles) model (for v s and uh 
considered here). 

• The non-coplanarity of the shock solution is also de- 
pendent upon the orientation of Bo, with the peak in B y 
decreasing as 6 — > 90° . 

• For pro-shock conditions considered here (tih = 
10 5 cm~ 3 , v s — 18km/s) the amount of dissipation inside 
the shock is relatively insensitive to the orientation of Bo, 
so in calculating molecular line emission or chemical models 
the shock problem is greatly simplified by suppressing B y . 
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